Extreme events in discrete nonlinear lattices 
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We perform statistical analysis on discrete nonlinear waves generated though modulational insta- 
bility in the context of the Salerno model that interpolates between the intergable Ablowitz-Ladik 
(AL) equation and the nonintegrable discrete nonlinear Schrodinger (DNLS) equation. We focus 
on extreme events in the form of discrete rogue or freak waves that may arise as a result of rapid 
coalescence of discrete breathers or other nonlinear interaction processes. We find power law depen- 
dence in the wave amplitude distribution accompanied by an enhanced probability for freak events 
close to the integrable limit of the equation. A characteristic peak in the extreme event probability 
appears that is attributed to the onset of interaction of the discrete solitons of the AL equation and 
the accompanied transition from the local to the global stochasticity monitored through the positive 
Lyapunov exponent of a nonlinear map. 

PACS numbers: 63.20.Ry; 47.20.Ky; 05.45+a 
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Introduction.- The motivation of the present work 
stems from observations of the sudden appearance of ex- 
tremely large amplitude sea waves referred to as rogue or 
freak waves [l[ . These waves appear very suddenly in rel- 
atively calm seas, reach amplitudes of over 20m and may 
destroy or sink small as well as large vessels [2|. The- 
oretical analysis of ocean freak waves has been linked 
to nonlinearities in the waver wave equations, studied 
though the nonlinear Schrodinger (NLS) equation and 
shown that the probability of their appearance is not in- 
significant Q. A scenario for freak wave generation in 
NLS is through a Benjamin- Feir (modulational) instabil- 
ity, resulting in self-focusing effects and subsequent for- 
mation of freak waves [4(. Modulational instability (MI) 
induces local exponential growth in the wave train ampli- 
tude [El, [f| that has been confirmed experimentally and 
numerically Q- 

Intriguingly, there are completely different physical 
systems that possess the required nonlinear characteris- 
tics which favour the appearance of rogue waves. Recent 
observation of optical rogue waves in a microstructured 
optical fiber was reported 0] in a regime near the thresh- 
old of soliton-fission supercontinuum generation, i.e., in 
a region where MI plays a key role in the dynamics. A 
generalized NLS equation was used successfully to model 
the generation of optical rogue waves while, additionally, 
control and manipulation of rogue soliton formation was 
also discussed [9[. The mechanism of the rogue waves 
creation, or, more generally of extreme events, has be- 
come an issue of principal interest in various other con- 
texts as well, since rogue waves can signal catastrophic 
phenomena such as an earthquake, a thunderstorm, or a 
severe financial crisis. Knowledge of the probability of 
occurrence of extreme events and the capability to pre- 
dict the time at which such an event may take place is of 
a great value. Such events are usually rare, and they ex- 



hibit " extreme- value" statistics, typically characterized 
by heavy-tailed probability distributions. Experimen- 
tal observation of optical rogue-wave-like fluctuations in 
fiber Raman amplifiers show that the probability distri- 
bution of their peak power follows a power law j 1 01 ] . 

In this work we focus on the discrete counterparts of 
rogue waves that may appear in nonlinear lattices as a 
result of discrete soliton or breather induction and their 
mutual interactions. Specifically we investigate the role 
of integrability in the formation of discrete rogue waves 
(DRW) and the resulting extreme event statistics. Their 
appearance may affect dramatically the physical systems. 
We use the Salerno model [III that through a unique 
parameter interpolates between a fully integrable discrete 
lattice, viz. the Ablowitz-Ladik (AL) lattice [l2j , and 
the nonintegrable DNLS equation [lj, Qjj]. One of the 
basic questions to be addressed below is the probability 
of occurrence of a DRW as a function of the degree of 
integrability of the lattice and thus study the role of the 
latter in the production of extreme lattice events [lo| . 

The Salerno model.- The Salerno model (SM) is given 
through the following set of equations 

i ~dT = ~^ + ^™| 2 )('^ n + 1 + ^n-l) - TlV'nlVn (1) 

where /i and 7 are two nonlinearity parameters. When 
(j, = the model becomes the DNLS equation while for 
7 = it reduces to AL. Several properties of the model 
such as inte grability [l6| and stability of localized trav- 
elling waves [171 Il8| | have been analyzed. Both the norm 
N and the Hamiltonian H of the model are conserved 
quantities. They are given by 
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It is also known that Eq. ([I]) exhibits MI, which may give 
rise to stationary, spatially localized solutions in the form 
of discrete breathers (DBs), i.e., periodic and spatially lo- 
calized nonlinear excitations [19j. The MI induced DBs 
appear in random lattice locations and may be mobile. 
High-amplitude DBs tend to absorb low-amplitude ones, 
resulting after some time in a small number of very high 
amplitude excitations, which may get pinned at a specific 
lattice site due to the Peierls-Nabarro potential barrier 
in nonintegrable lattices 20]. In general, high-amplitude 
DBs are virtual bottlenecks which slow down the relax- 
ation processes in nonlinear lattices [2l|, , and it has 
been proposed that they may serve as models for freak 
waves [23| . The development of MI in Eq. (p} can be an- 
alyzed with the linear stability analysis of its the plane 
wave solutions perturbed by small phase and amplitude 
perturbations [2J] . The interplay between the on-site and 
intersite nonlinear terms (i.e., according to the variation 
of their relative strength through [i and 7), may change 
MI properties and, consequently, the conditions for the 
DBs to exist in the lattice 25]. The SM has recently 
found applications in modelling Bose-Einstein conden- 
sates of dipolar atoms in a strong periodic potential [26| , 
dilute Bose-Einstein condensates trapped in a periodic 
potential [27| , and even biological systems [28[ . 

For later convenience in the numerical simulation, the 
variables ip n in Eq. ([T]) are rescaled as ip n — £, n / \JJi, so 
that in terms of the dynamic equations read 

= _ (1 + |£ n(t) |2 )( £ n+i + u _ x) _ n^)i 2 £ n ,(4) 



df 



where T — 7///. Therefore, the whole two-dimensional 
parameter space (7, fx) can be scaled by // = 1, leaving 
7 as a free parameter. With that scaling we may go as 
close to the DNLS limit as we want to, by simply let T to 
attain very large values. However, the exact DNLS limit 
li = has to be calculated separately. 

Statistics of extreme events.- We integrate numerically 
the system of Eqs. @ with periodic boundary conditions 
using a sixth order Runge-Kutta algorithm with fixed 
time-stepping At = 10~ 4 . We started simulations with 
different initial conditions (the plane wave, uniform back- 
ground with white noise and Gaussian noise) which gave 
similar results. Here we present calculations in which 
the initial condition is uniform, = 1 for any n, with 
the addition of a small amount of white noise to accel- 
erate the development of the MI. The uniform solution 
is chosen in the interval where it is known from linear 
stability analysis that it is unstable. By varying the non- 
linearity parameters we identify broadly three regimes of 
DRWs that are shown as spatiotemporal evolutions in 
Fig. Q] For the purely integrable AL lattice (r = 0) DBs 
are mobile and essentially noninteracting; as a result we 
do not observe significant formation of high DRWs (Fig. 
QJi). In the vicinity of the AL limit, i.e. for small T 
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FIG. 1: Evolution of the scaled amplitudes for a lattice 
of size N = 101, with V (/i and 7 in the DNLS case), is shown 
on the figure. The initial conditions for all cases are f n = 1 
for any n (uniform) plus a small amount of white noise. 



(r ~ 0.1), there is an onset of weak interaction of the 
localized modes of the AL lattice leading to a significant 
increase in DRW formation that are mobile (Fig. [T}),c). 
In this regime the DBs are highly mobile indicating that 
DB merging could be responsible for creation of high- 
amplitude localized waves. For T >> 0.1 on the other 
hand, DNLS-type behavior dominates the SM and local- 
ized structures that are initially created through the MI 
become easily trapped in the lattice (Fig. QJ1). 

The three regimes mentioned previously are probed 
by calculating the time-averaged height distributions Ph- 
We first define the forward (backward) height at the 
n— th site as the difference between two successive min- 
imum (maximum) and maximum (minimum) values of 
|£„(t)|. We use then both the forward and the backward 
heights for the calculation of the local height distribution; 
after spatial averaging the latter results in the height 
probability densities (HPDs) shown in Fig. 2. We note 
that the tails of the HPDs are related to extreme events 
and the appearance of DRWs. For T finite the HPDs are 
sharply peaked but have extended tails indicating that 
extreme events are more than several times as large as 
the mean distribution height. In the DNLS limit (r ^> 1) 
the obtained HPD is very close to the Rayleigh distribu- 
tion whose tails decay very fast [2^] , indicating negligible 
probability for the occurrence of extreme events (dotted 
curve in Fig. [2]). In all the other cases the decay of the 
tails of the HPDs is much slower. 

In order to probe further the onset of extreme dis- 
crete events we employ the practice used in water waves 
and define a DRW as one that has a height greater than 
h t h = 2.2/i s , with h s being the significant wave height. 
The latter is defined as the average height of the one- 
third higher waves in the height distribution. As a result, 
the probability of occurrence of extreme DRW events 
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FIG. 2: The normalized height probability density Ph(h) for 
several values of F and for the DNLS limit (with 7 = 6). 
The line with slope —1 is added to assist comparisons and 
corresponds to Ph ~ 1/h. Approximately vertical drop cor- 
responds to the DNLS limit with an exponential tail. The 
increase of F(^ = 1) leads to the decrease of the slope and 
appearance of plateau on the Ph curve; the latter increases 
the extreme event probability leading a maximum at T = 0.07 
(Fig. 3.) 

Pee = Ph(h > hth) is obtained by integration of the 
(normalized) HPD from h = hth up to infinity. By eval- 
uating several HPDs as those in Fig. [2] we may estimate 
the probability of occurrence of DRWs P ee as a function 
of the parameter T (the results are shown in Fig. [3J . We 
note that the probability for the occurrence of a DRW 
has a certain value in the AL case, subsequently peaks for 
small values of T and decays precipitously when T » 1. 
This behavior of the probability P ee is compatible with 
the DB picture outlined earlier, viz. in the very weakly 
nonintegrable regime the AL modes may interact leading 
to DB fusion and DRW generation. On the other hand, 
as nonintegrability becomes stronger, the scattering of 
the AL modes is more chaotic leading to a suppression 
of DRW formation. 

Map approach.- In order to probe deeper on the for- 
mation of DRWs we substitute ip n = cj) n exp(— iojt) into 
Eq. (pj, with <j) n a real- valued function of the lattice site 
n, and obtain the stationary equation 

LU(j> n + (1 + fl\<P n \ 2 )((j> n+1 + (t> n - X ) + 7 |0„| 2 0n = 0, (5) 

which can be transformed in the two-dimensional map 

x n +i = -— ^-x n - y n , y n+1 = x n , (b) 

where we have defined x n = 4> n and y n = <fi n -i- Eqs . 
([6|) represent a real analytic area-preserving map [LsL I30I ] 
with the lattice index n playing the role of discrete 'time'. 

The phase portraits of the map Eq. for several 
T-values are shown in Fig. |4j In the AL limit, the 
phase space consists of perfectly disconnected separatri- 
ces while for non-zero T, the stable and unstable man- 
ifolds intersect transversely, resulting in the generation 



FIG. 3: The normalized probability P ee — Ph(h > hth) f° r 
the occurrence of extreme events as a function of the integra- 
bility parameter F. All data present averaged results of five 
numerical measurements differing in the initial conditions. 



of a homoclinic tangle. With increasing T the motion 
near separatrices becomes exceedingly complicated and 
the trajectories wander irregularly before approaching an 
attracting set (Figs. 4b and 4c). Moreover, for any r ^ 0, 
the position of separatrices in phase space changes in 
time, resulting in overlapping of neighboring separatrices 
and diffusion in those regions which have been traversed 
by a separatrix. The sharp peak of the probability of 
occurrence of extreme events P ee (h > hth) in the SM 
(Fig. [3J can be associated with the opening of a stochas- 
ticity web, when orbits fast explore all extended narrow 
stochasticity regions leading to an anomalous relaxation 
phase [lj| [29[ • This event signs the transition from the lo- 
cal to global stochasticity [3i( in SM. On the other hand, 
the decrease of P ee (h > hth) for larger T's is related to 
the increasingly longer trapping time in more developed 
stochasticity region. 

The Melnikov analysis in the SM [l8| shows that the 
magnitude of the separatrix splitting and the consequent 
development of stochasticity depends on the T/\ui\ ra- 
tio. The conjecture that P ee is associated with the com- 
plexity of the phase portraits of the corresponding maps 
implies that P ee should also depend on the r/|w| ra- 
tio. In our case |u;| is related to the modulation fre- 
quency of the initially uniform solution U with the rela- 
tion \oj\ = (7 + 2/i)[/ 2 + 2, which, through the MI process 
it transformed into a train of localized DB-like config- 
urations. We have checked numerically that for fixed 
ratio T/\uj\ and different values of U and T we obtain the 
same HPD. As a consequence, the probability of extreme 
events P ee as a function of the r/|w| is qualitatively the 
same with that of P ee as a function of T shown in Fig. [3] 

The degree of nonintegrability in the SM model can 
be quantified by calculating the Lyapunov exponents of 
the corresponding maps [32|]. We have thus calculated 
the maximum Lyapunov exponent L [3l| for the map 
Eq. ©, for the parameters used in the calculation of 
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FIG. 4: Orbits started at different initial positions in the 
neighborhood of map origin and corresponding the one- 
dimensional Liapunov exponents. 

the phase portrait shown in the left panels of Fig. [4j ft 
is observed that homoclinic orbits which correspond to 
perfect separatrices are characterized by vanishing Lya- 
punov exponent (Fig. 2^). With increasing stochasticity, 
L tends to a finite positive value which generally depends 
on the values of the parameters and the initial conditions 



(Figs. Ha and |4p). 

Conclusions.- The probability of occurrence of extreme 
events P ee in the SM results from the competition be- 
tween the self-focusing and the energy transport mecha- 
nisms which are implicitly correlated with the degree of 
integrability of the model [IB]. Through modulational 
instability and starting from a slightly perturbed uni- 
form background we can generate high-amplitude local- 
ized moving structures of the DB type that lead to the 
formation of extreme events of DRW type. Depending 
on their number, amplitude and life-time, they may pre- 
vent of facilitate the energy flow in the lattice, affect- 
ing thus the probability of extreme event formation P ee . 
We find that the latter probability depends strongly on 
r that affects the degree of integrability of the lattice: 
DRW are much more probable very close to the inte- 
grable SM limit rather than in the nonintegrable one. We 
find a resonance-like maximum in P ee (T) that, through 
a nonlinear map approach, is linked to separatrix break- 
ing and the onset of global stochasticity. This regime 
corresponds physically to weak interaction between the 
quasi-integrable modes of the system. 
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